
from os import listdir
from math import sqrt, ceil, floor, pow, pi
from matplotlib import pyplot

def Lboundary(x, lx):
    res=x
    while x>=lx:
        x-=lx
    while x<0:
        x+=lx
    return x

def boundary(pos, lx,ly,lz):
    pos[0]=Lboundary(pos[0],lx)
    pos[1]=Lboundary(pos[1],ly)
    #pos[2]=Lboundary(pos[2],lz)
    return pos

def gr(pos, side, numbins):
    
    maxbin = numbins # number of bins
    sideh = side/2.0
    dr = float(sideh/maxbin) # bin width
    hist = [0]*(maxbin+1) # initialize list
    count = 0
            
    # read atom coordinates from PDB
    atoms = []
    for coords in pos:   
        atoms.append((coords[0], coords[1], coords[2]))             
        
    # loop over particle pairs
    npart = len(atoms)
    
    for i in range(npart):
        
        xi = (atoms[i])[0]
        yi = (atoms[i])[1]
        zi = (atoms[i])[2]
        
        for j in range(i+1, npart):
            xx = xi - (atoms[j])[0]
            yy = yi - (atoms[j])[1]
            zz = zi - (atoms[j])[2]
            
            # minimum image
            if (xx < -sideh):   xx = xx + side
            if (xx > sideh):    xx = xx - side
            if (yy < -sideh):   yy = yy + side
            if (yy > sideh):    yy = yy - side
            if (zz < -sideh):   zz = zz + side
            if (zz > sideh):    zz = zz - side
            
            # distance between i and j
            rd  = xx * xx + yy * yy + zz * zz
            rij = sqrt(rd)
            
            bin = int(ceil(rij/dr)) # determine in which bin the distance falls
            if (bin <= maxbin):
                hist[bin] += 1
    
    # normalize
    phi = npart/pow(side, 3.0) # number density (N*V)
    norm = 2.0 * pi * dr * phi  * npart
    
    x=[]
    y=[]
    for i in range(1, maxbin+1):
        rrr = (i - 0.5) * dr
        val = hist[i]/ norm / ((rrr * rrr) + (dr * dr) / 12.0)
        x.append(rrr)
        y.append(val)
    
    #return rdf
    return [x,y]

def grXY(pos, side, numbins):
    
    maxbin = numbins # number of bins
    sideh = side/2.0
    dr = float(sideh/maxbin) # bin width
    hist = [0]*(maxbin+1) # initialize list
    count = 0
            
    # read atom coordinates from PDB
    atoms = []
    for coords in pos:   
        atoms.append((coords[0], coords[1], coords[2]))             
        
    # loop over particle pairs
    npart = len(atoms)
    
    for i in range(npart):
        
        xi = (atoms[i])[0]
        yi = (atoms[i])[1]
        #zi = (atoms[i])[2]
        
        for j in range(i+1, npart):
            xx = xi - (atoms[j])[0]
            yy = yi - (atoms[j])[1]
            #zz = zi - (atoms[j])[2]
            
            # minimum image
            if (xx < -sideh):   xx = xx + side
            if (xx > sideh):    xx = xx - side
            if (yy < -sideh):   yy = yy + side
            if (yy > sideh):    yy = yy - side
            #if (zz < -sideh):   zz = zz + side
            #if (zz > sideh):    zz = zz - side
            
            # distance between i and j
            #rd  = xx * xx + yy * yy + zz * zz
            rd  = xx * xx + yy * yy
            rij = sqrt(rd)
            
            bin = int(ceil(rij/dr)) # determine in which bin the distance falls
            if (bin <= maxbin):
                hist[bin] += 1
    
    # normalize
    phi = npart/pow(side, 3.0) # number density (N*V)
    norm = 2.0 * pi * dr * phi  * npart
    
    x=[]
    y=[]
    for i in range(1, maxbin+1):
        rrr = (i - 0.5) * dr
        #val = hist[i]/ norm / ((rrr * rrr) + (dr * dr) / 12.0)
        val = hist[i]/ norm / ((rrr ) + (dr ) / 12.0)/(0.5*side)
        x.append(rrr)
        y.append(val)
    
    #return rdf
    return [x,y]

def importConfig(name):
    f = open(name,'r')
    tx = f.readline()
    tx = tx.replace('\n','')
    val = tx.split(' ')
    aux=[]
    for v in val:
        if v!='':
            aux.append(v)
    t=float(aux[1])
    tcorr=float(aux[3])
    icorr=int(aux[5])
    npart=int(aux[7])
    lx=float(aux[9])
    ly=float(aux[10])
    lz=float(aux[11])
    tx = f.readline()
    tx = tx.replace('\n','')
    val = tx.split(' ')
    aux=[]
    for v in val:
        if v!='':
            aux.append(v)
    sigma=[float(aux[1]),float(aux[2]),float(aux[3]),float(aux[4]),float(aux[5]),float(aux[6])]
    pos=[]
    vel=[]
    force=[]
    diameter=[]
    
    for i in range(npart):
        tx = f.readline()
        tx = tx.replace('\n','')
        val = tx.split(' ')
        aux=[]
        for v in val:
            if v!='':
                aux.append(v)
    
        xpos=boundary([float(aux[0]),float(aux[1]),float(aux[2])], lx,ly,lz)
        pos.append(xpos)
        vel.append([float(aux[3]),float(aux[4]),float(aux[5])])
        force.append([float(aux[6]),float(aux[7]),float(aux[8])])
        diameter.append(float(aux[9]))
    return [pos,vel,force,diameter]

name = '/home/jonk/eclipse/analisis/salida.dat'
[pos,vel,force,diameter] = importConfig(name)
side=27.395611
numbins = 100
[x,y] = gr(pos, side, numbins)
pyplot.plot(x,y,'r-')
[x,y] = grXY(pos, side, 2*numbins)
pyplot.plot(x,y,'b-')
pyplot.show()

